josef — Aug 5, 2014, 11:04 PM
## # Computational Exercise 3.1: Estimating and plotting a selection surface
## ### 2012 Stevan J. Arnold
## ## Get your data into R using the following steps.
## Save the file radix5.txt to your desktop or folder of choice
## In the following example, we have created a folder called ?R wd? and placed the data file in it
## Now, set your working directory using a statement comparable to this one
## Be sure and use the / not the \ spacing convention!
#setwd('set working directory here')
## Now read in your data
thamnophis = read.table("radix5.txt",header=T)
#Print out the data frame to make sure it looks OK
head(thamnophis)
LITTER NAME SPEED BODY TAIL
1 761 1 0.16207 154 77
2 761 3 0.09019 149 70
3 761 4 0.10823 153 69
4 761 5 0.27609 154 73
5 761 6 0.13532 150 72
6 761 8 0.20325 149 68
## ## Stage 1: Estimating the selection gradients
## Use the following statements to create standardized versions of BODY, TAIL, and TIME
new.body=mean(thamnophis$BODY)-thamnophis$BODY
new.tail=mean(thamnophis$TAIL)-thamnophis$TAIL
new.speed=thamnophis$SPEED/sd(thamnophis$SPEED)
## What do these transformations do to the mean and standard deviation of the variables?
## Find out by using statements like 'mean(new.body), sd(new.body)'.
## Let's begin by fitting a plane to the transformed data
model1 <-lm(new.speed~new.body + new.tail)
## print out the coefficients and statistics of the fit
## this output will give us our best estimates of the directional selection gradients for
## new.body and new.tail
summary(model1)
Call:
lm(formula = new.speed ~ new.body + new.tail)
Residuals:
Min 1Q Median 3Q Max
-2.2202 -0.6755 0.0165 0.7030 2.4108
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.62196 0.08393 31.24 <2e-16 ***
new.body -0.02435 0.02497 -0.98 0.33
new.tail -0.00131 0.02800 -0.05 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1 on 140 degrees of freedom
Multiple R-squared: 0.00684, Adjusted R-squared: -0.00734
F-statistic: 0.482 on 2 and 140 DF, p-value: 0.618
## Now, fit a full quadratic model to the data, remembering to use the factor of 0.5 for the
## stabilizing selection gradients
## Create a new variable called prod which is the product of new.body and new.tail
prod = new.body*new.tail
## Then fit the quadratic model
model2 <- lm(new.speed ~ new.body + new.tail + I(0.5*new.body^2) + I(0.5*new.tail^2) + prod )
summary(model2)
Call:
lm(formula = new.speed ~ new.body + new.tail + I(0.5 * new.body^2) +
I(0.5 * new.tail^2) + prod)
Residuals:
Min 1Q Median 3Q Max
-2.2180 -0.6022 0.0418 0.6727 2.4024
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.62844 0.11395 23.07 <2e-16 ***
new.body -0.02759 0.02489 -1.11 0.270
new.tail 0.00809 0.02800 0.29 0.773
I(0.5 * new.body^2) -0.00243 0.00914 -0.27 0.791
I(0.5 * new.tail^2) -0.00171 0.01336 -0.13 0.899
prod 0.02031 0.00789 2.57 0.011 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.99 on 137 degrees of freedom
Multiple R-squared: 0.0537, Adjusted R-squared: 0.0192
F-statistic: 1.56 on 5 and 137 DF, p-value: 0.177
## In summary, using z1 for body and z2 for tail,
## we have from the results above, the following values for selection gradients
beta1= coef(model1)[2]
beta2= coef(model1)[3]
gamma11 = coef(model2)[4]
gamma22 = coef(model2)[5]
gamma12 = coef(model2)[6]
## ## Stage 2: Plotting the selection surface, an ISS
## Set the parameters for a full quadratic fitness surface for two traits
z1 = new.body
z2 = new.tail
## Set the value for the intercept of the surface when z1=z2=0
alpha=1
## Define a function called fit that will compute the value of fitness as a function of z1 and z2
fit <- function(z1,z2, alpha, beta1, beta2, gamma11, gamma22, gamma12) alpha + (beta1*z1) + (beta2*z2) + (gamma11*0.5*(z1^2)) + (gamma22*0.5*(z2^2)) + (gamma12*(z1 * z2))
## Define a series of values for z1 and z2 that will be used to compute the value of relative fitness on the surface
x <- seq(-2, 2, length = 30)
y <- seq(-2, 2, length = 30)
## Compute the surface values of relative fitness using the x-y grid of values for z1 and
## z2 for later use by the surface plotting function called persp
## We could use the fit function to compute values of fitness but instead we will use a #function, called outer, that is more compatible with persp
## The function outer has some specific requirements so we oblige by writing our fitness
## function in the following form
z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( -0.024351*a) + (-0.001307*b) + (-0.002425*0.5*(a^2)) + (-0.001705*0.5*(b^2)) + (0.020315 *(a * b)))
## Define two variables that give the number of rows and cols in z
nrz <- nrow(z)
ncz <- ncol(z)
## Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("yellow", "orange") )
## Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
## Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
## Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
## Finally, plot the ISS
par(bg = "white")
persp(x, y, z, col=color[facetcol], phi=30, theta=-30)
## For another view of the surface, do a simple contour plot
par(pty ="s")
contour(x, y, z)
## Let's plot the eigenvectors of the gamma-matrix on this contour plot
## First, let?s define the matrix
gamma = matrix(c(gamma11, gamma12, gamma12, gamma22), nrow=2, ncol=2)
gamma
[,1] [,2]
[1,] -0.002425 0.020315
[2,] 0.020315 -0.001705
eigen(gamma)
$values
[1] 0.01825 -0.02238
$vectors
[,1] [,2]
[1,] -0.7008 -0.7133
[2,] -0.7133 0.7008
## Next, let?s setup the beta vector and check to see if its values are correct
beta = matrix(c(beta1, beta2), nrow=2, ncol=1)
beta
[,1]
[1,] -0.024351
[2,] -0.001307
## Using gamma and beta, we can solve for the stationary point on the surface using
## expression(11) from Phillips & Arnold 1989
z.zero = -(solve(gamma)) %*% beta
## Plot the stationary point on the surface
points(z.zero[1,1], z.zero[2,1])
## Finally, plot the eigenvectors on the surface
## First one eigenvector
delx1 = z.zero[1] + 6*eigen(gamma)$vectors[1,2]
delx2 = z.zero[2] + 6*eigen(gamma)$vectors[2,2]
delx11 = z.zero[1] - 6*eigen(gamma)$vectors[1,2]
delx22 = z.zero[2] - 6*eigen(gamma)$vectors[2,2]
x.values=c(delx1, z.zero[1], delx11)
y.values=c(delx2, z.zero[2], delx22)
lines(x.values, y.values, lty=2,lwd=2)
## then the other
delx1 = z.zero[1] + 6*eigen(gamma)$vectors[1,1]
delx2 = z.zero[2] + 6*eigen(gamma)$vectors[2,1]
delx11 = z.zero[1] - 6*eigen(gamma)$vectors[1,1]
delx22 = z.zero[2] - 6*eigen(gamma)$vectors[2,1]
x.values=c (delx1, z.zero[1], delx11)
y.values=c(delx2, z.zero[2], delx22)
lines(x.values, y.values, lty=2, lwd=2)
## Which of these eigenvectors is gamma max?
## Now, let?s estimate the omega-matrix, approximating it as the negative inverse of the
## gamma-matrix
omega=-(solve(gamma))
omega
[,1] [,2]
[1,] -4.173 -49.723
[2,] -49.723 -5.936
## It's eigenvectors should be the same as those of the gamma-matrix
## Let's check
eigen(omega)
$values
[1] 44.68 -54.79
$vectors
[,1] [,2]
[1,] -0.7133 0.7008
[2,] 0.7008 0.7133
## What is the interpretation of these eigenvalues on the contour plot?
## ## Stage 3: Can you build a loop that will show the contour plot while bootstrapping #over the snake sample? Here?s an example of script that does that!
## We begin by writing a function that we will use during bootstrapping to estimate the coefficients that describe the surface for a particular boot sample
est.coeff=function(y,x1,x2){
prod=x1*x2
m1=lm(y~x1+x2)
m2=lm(y~x1+x2+I(0.5*x1^2)+I(0.5*x2^2)+prod)
c(m1$coeff[2],m1$coeff[3],m2$coeff[4],m2$coeff[5],m2$coeff[6])
}
## Next, we program a perspective plot function
x <- seq(-2, 2, length = 30)
y <- seq(-2, 2, length = 30)
## Compute the surface values of relative fitness using the x-y grid of values for z1 and
## z2 for later use by the surface plotting function called persp
## We could use the fit function to compute values of fitness but instead we will use a #function, called outer, that is more compatible with persp
## The function outer has some specific requirements so we oblige by writing our fitness
## function in the following form
## Write a function that plots the perspective plot for each bootstrap
persp.boot=function(b1,b2,y1,y2,y12,x=seq(-10,10,length=30),y=seq(-10,10,length=30)){
z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( b1*a) + (b2*b) + (y1*0.5*(a^2)) + (y2*0.5*(b^2)) + (y12*(a * b)))
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette( c("yellow", "orange") )
nbcol <- 100
color <- jet.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)
par(bg = "white")
persp(x, y, z, col=color[facetcol],phi=30,theta=-30)
}
## Test it out
persp.boot(-2,2,1,1,1)
## Bootstrap over individuals
par(ask=FALSE)
n=length(new.body)
for (i in 1:100){
samp=sample(1:n,n,replace=TRUE)
boot.speed=new.speed[samp]
boot.body=new.body[samp]
boot.tail=new.tail[samp]
boot.coeff=est.coeff(boot.speed,boot.body,boot.tail)
persp.boot(boot.coeff[1],boot.coeff[2],boot.coeff[3],boot.coeff[4],boot.coeff[5],x=seq(-10,10,length=10),y=seq(-10,10,length=10))
Sys.sleep(0.1)
}
## ## Stage 4: Let's build a loop that will show the contour plot while bootstrapping over the snake sample while rotating the surface about this vertical axis?
## First, let's make the perspective plot rotate
k=100
for (i in 1:k){
persp(x, y, z, col=color[facetcol],phi=30,theta=-(30 +(3.3*i)))
}
## To make the plot rotate while bootstrapping we first
## write a new function that plots the perspective plot for each bootstrap
i=1
persp.boot=function(b1,b2,y1,y2,y12,x=seq(-10,10,length=30),y=seq(-10,10,length=30)){
z <- outer(x, y, function(a, b, alpha, beta1, beta2, gamma11, gamma22, gamma12) 1 + ( b1*a) + (b2*b) + (y1*0.5*(a^2)) + (y2*0.5*(b^2)) + (y12*(a * b)))
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette( c("yellow", "orange") )
nbcol <- 100
color <- jet.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)
par(bg = "white")
persp(x, y, z, col=color[facetcol],phi=30,theta=-30+(3.3*i))
}
## Not let's run the bootstrapping during the rotation
par(ask=FALSE)
n = length(new.body)
k=100
for (i in 1:k){
samp=sample(1:n,n,replace=TRUE)
boot.speed=new.speed[samp]
boot.body=new.body[samp]
boot.tail=new.tail[samp]
boot.coeff=est.coeff(boot.speed,boot.body,boot.tail)
persp.boot(boot.coeff[1],boot.coeff[2],boot.coeff[3],boot.coeff[4],boot.coeff[5],x=seq(-10,10,length=10),y=seq(-10,10,length=10))
Sys.sleep(0.2)
}